ohi logo
OHI-Northeast | OHI Science | Citation policy


#Data Source

Downloaded: December 14, 2018 (emailed to us by Jeffrey Vieser at NMFS)

Description: Records of Bmsy and Fmsy estimates from stock assessments conducted in the greater Northeast region

Time range: 2004 - 2018

Format: Tabular


1 Data wrangling

Load raw data

Clean data

Save data

Span plot

What species have multiple stocks? These ones will give us a tough time when matching up catch to the stocks.

## # A tibble: 20 x 3
##    stock                             species        location               
##    <chr>                             <chr>          <chr>                  
##  1 Atlantic cod - Georges Bank       Atlantic cod   Georges Bank           
##  2 Atlantic cod - Gulf of Maine      Atlantic cod   Gulf of Maine          
##  3 Haddock - Georges Bank            Haddock        Georges Bank           
##  4 Haddock - Gulf of Maine           Haddock        Gulf of Maine          
##  5 Windowpane - Gulf of Maine / Geo… Windowpane     Gulf of Maine / George…
##  6 Windowpane - Southern New Englan… Windowpane     Southern New England /…
##  7 Winter flounder - Gulf of Maine   Winter flound… Gulf of Maine          
##  8 Winter flounder - Southern New E… Winter flound… Southern New England /…
##  9 Yellowtail flounder - Cape Cod /… Yellowtail fl… Cape Cod / Gulf of Mai…
## 10 Yellowtail flounder - Georges Ba… Yellowtail fl… Georges Bank           
## 11 Yellowtail flounder - Southern N… Yellowtail fl… Southern New England /…
## 12 Winter flounder - Georges Bank    Winter flound… Georges Bank           
## 13 Red hake - Gulf of Maine / North… Red hake       Gulf of Maine / Northe…
## 14 Red hake - Southern Georges Bank… Red hake       Southern Georges Bank …
## 15 Silver hake - Gulf of Maine / No… Silver hake    Gulf of Maine / Northe…
## 16 Silver hake - Southern Georges B… Silver hake    Southern Georges Bank …
## 17 Atlantic cod - Eastern Georges B… Atlantic cod   Eastern Georges Bank   
## 18 Haddock - Eastern Georges Bank    Haddock        Eastern Georges Bank   
## 19 Goosefish - Gulf of Maine / Nort… Goosefish      Gulf of Maine / Northe…
## 20 Goosefish - Southern Georges Ban… Goosefish      Southern Georges Bank …

Looks like there are just 8.

2 Gapfill

Use a linear model to fill years between assessments.

2.1 Calculate stock scores

Convert these value to stock scores.

#grab just B/Bmsy data
  nmfs_b_bmsy <- data_gf %>%
    filter(indicator == "b_bmsy") %>%
    select(stock, year, value)

#grab just F/Fmsy data  
  nmfs_f_fmsy <- data_gf %>%
    filter(indicator == "f_fmsy") %>%
    select(stock, year, value)
  ### 
  b_bmsy_underexploit_penalty <- 0.25
  b_bmsy_underexploit_thresh  <- 3.00
  f_fmsy_underfishing_penalty <- 0.25
  f_fmsy_overfishing_thresh   <- 2.00
  
  ### Apply rolling mean to F/Fmsy
  ## Why do we do this? because B is a less sensitive metric (relies of biological processes) and F can fluctuate pretty easily because it is really just a mgmt decision.
  nmfs_f_fmsy <- nmfs_f_fmsy %>%
    arrange(stock, year) %>%
    group_by(stock) %>%
    mutate(f_fmsy_rollmean = zoo::rollmean(value, k = 3, align = 'right', fill = NA)) %>%
    ungroup() %>%
    select(-value) %>%
    rename(value = f_fmsy_rollmean)
  
  stock_status_layers <- nmfs_b_bmsy %>%
    full_join(nmfs_f_fmsy) %>%
    spread(indicator, value) 
  
########################################################.
##### run each fishery through the Kobe plot calcs #####
########################################################.
### * ram_b_bmsy, ram_f_fmsy
  
  
### Function for converting B/Bmsy values into a 0 - 1 score
  rescale_bprime_crit <- function(fish_stat_df,
                                  bmax, bmax_val) {
    
    ###using NOAA's limits here
    overfished_th  <- 0.8
    ### 
    underfished_th <- 1.2
    
    bmax_adj <- (bmax - underfished_th) / (1 - bmax_val) + underfished_th
    ### this is used to create a 'virtual' B/Bmsy max where score drops
    ### to zero.  If bmax_val == 0, this is bmax; if bmax_val > 0, bmax_adj
    ### extends beyond bmax, to create a gradient where bmax_val occurs at bmax.
    
    fish_stat_df <- fish_stat_df %>%
      # group_by(stock) %>% ### grouping by stock will set b_max by max per stock, instead of max overall
      mutate(b_max     = max(b_bmsy, na.rm = TRUE)) %>%
      ungroup() %>%
      mutate(bPrime = NA,
             bPrime = ifelse(b_bmsy < overfished_th,  ### overfished stock
                             b_bmsy / overfished_th,
                             bPrime),
             bPrime = ifelse(b_bmsy >= overfished_th & b_bmsy < underfished_th,
                             1,                       ### appropriately fished stock
                             bPrime),
             bPrime = ifelse(b_bmsy >= underfished_th,
                             (bmax_adj - b_bmsy) / (bmax_adj - underfished_th), ### underfished stock
                             bPrime),
             bPrime = ifelse(bPrime < 0, 0, bPrime))
    return(fish_stat_df)
  }
  
  
  ### Function to create vertical gradient based on distance from
  ### ideal F/Fmsy value to actual F/Fmsy value
  f_gradient <- function(f, over_f, under_f, fmax, fmin_val) {
    x <- ifelse(f < over_f & f > under_f, 1, NA)
    x <- ifelse(f <= under_f, (f * (1 - fmin_val) / under_f + fmin_val), x)
    x <- ifelse(f >= over_f,  (fmax - f) / (fmax - over_f), x)
    x <- ifelse(f > fmax, NA, x)
    return(x)
  }
  
  ### Function to convert F/Fmsy values into 0 - 1 score
  rescale_fprime_crit <- function(fish_stat_df,
                                  fmax, fmin_val) {
    
    ### params - taken from BC but changed Bcrit to 0 instead of 0.4:
    Bcrit <- 0.5; overfished_th <- 0.8
    ### underfishing_th is set to the idea "1/3 for the birds":
    underfishing_th <- 0.66; overfishing_th  <- 1.2
    
    bcritslope = 1 / (overfished_th - Bcrit)
    ### connecting from (Bcrit, 0) to (overfished_th, 1)
    
    fish_stat_df <- fish_stat_df %>%
      mutate(fPrime = ifelse(b_bmsy < overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy + (overfished_th - b_bmsy) * bcritslope,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             NA),
             fPrime = ifelse(b_bmsy >= overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             fPrime),
             fPrime = ifelse(is.na(fPrime), 0, fPrime), ### fill zeros everywhere unscored
             fPrime = ifelse(is.na(f_fmsy), NA, fPrime) ### but if no f_fmsy, reset to NA
      )
    return(fish_stat_df)
  }
  
  stock_status_df <- stock_status_layers %>%
    rescale_bprime_crit(bmax     = b_bmsy_underexploit_thresh,
                        bmax_val = b_bmsy_underexploit_penalty) %>%
    rescale_fprime_crit(fmax     = f_fmsy_overfishing_thresh,
                        fmin_val = f_fmsy_underfishing_penalty) %>%
    mutate(x_prod = ifelse(!is.na(fPrime), (fPrime * bPrime), bPrime),
           basis  = case_when(
             !is.na(fPrime) & !is.na(bPrime) ~ 'F_Fmsy, B_Bmsy',
             is.na(fPrime)  & !is.na(bPrime) ~ 'B_Bmsy only',
             is.na(bPrime)  & !is.na(fPrime) ~ 'F_Fmsy only'
           )) %>%
    dplyr::select(year, stock,
                  score = x_prod,
                  basis,
                  bPrime, fPrime,
                  b_bmsy, f_fmsy) 

Take the average scores for the 5 species with sub-stocks. Just doing a group_by with species and averaging scores will do this for all species (but only change the values for these sub-species).

Remove stocks with just F_Fmsy - these stocks have no stock scores since we don’t have a way to get a score from just F/Fmsy. Also roll all values forward to 2017.

2.2 Visualize

kobe_plot <- function(sp){
  
 ss_df <- stock_scores %>%
    filter(stock == sp) %>%
    arrange(year) %>%
    mutate(last_bbmsy = last(b_bmsy),
           last_ffmsy = last(f_fmsy),
           last_datayear = last(year)) %>%
   ungroup()

generate_kobe_df <- function(f_fmsy_max = 2.5,
                             b_bmsy_max = 3.0,
                             reso       = 0.01,
                             bmax_val   = 0,
                             fmin_val   = 0,
                             weighting_b = 1) {

  kobe_raw <- data.frame(stock  = 1,
                     f_fmsy = rep(seq(0, f_fmsy_max, reso), each  = round(b_bmsy_max/reso) + 1),
                     b_bmsy = rep(seq(0, b_bmsy_max, reso), times = round(f_fmsy_max/reso) + 1))

  kobe <- kobe_raw %>%
    rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
                        bmax_val = bmax_val) %>%
    rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
                        fmin_val = fmin_val) %>%
    mutate(x_geom  = (fPrime * bPrime),
           x_arith = (fPrime + bPrime) / 2)

  return(kobe)
}

bbmsy_lim <- max(round(max(ss_df$b_bmsy, na.rm = TRUE) + .1, 1), 3)
ffmsy_lim <- max(round(max(ss_df$f_fmsy, na.rm = TRUE) + .1, 1), 2.5)
  
kobe_df <- generate_kobe_df(f_fmsy_max = ffmsy_lim,
                           b_bmsy_max = bbmsy_lim,
                           bmax_val = .25,
                           fmin_val = .25)


kobe_stock_plot <- ggplot(data = kobe_df, aes(x = b_bmsy, y = f_fmsy)) +
    theme_bw() +
    geom_raster(alpha = .8, aes(fill = x_geom)) +
    scale_fill_distiller(palette = 'RdYlGn', direction = 1) +
    labs(title = as.character(sp),
         x = 'B/Bmsy',
         y = 'F/Fmsy',
         fill = "Stock score") +
    geom_path(data = ss_df, 
              show.legend = FALSE,
              aes(x = b_bmsy, y = f_fmsy, group = sp),
              color = 'grey30') +
    geom_point(data = ss_df, 
               show.legend = FALSE,
              aes(x = last_bbmsy, y = last_ffmsy)) +
    geom_text(data = ss_df %>%
                mutate(year = ifelse(year/5 == round(year/5) | year == last_datayear, year, NA)), 
              aes(x = b_bmsy, y = f_fmsy, label = year), 
              hjust = 0, nudge_x = .05, size = 2)

return(kobe_stock_plot)
}

Apply function to all species in the stock_scores data frame.